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Energy level statistics of mixed quantum systems whose classical dynamics is partly 
regular and partly chaotic have been intensively studied over the past decade (see 
[!] and references therein), and this subject is still much less theoretically under- 
stood than the level statistics of the two extreme cases, namely completely chaotic 
(hyperbolic) systems 0, |3|, and integrable systems ||. However, it is believed that 
mixed systems, for example hydrogen atom in strong magnetic field [||], are generic 
in nature, at least among dynamical systems with few degrees of freedom. Although 
Berry and Robnik have developed a semiclassical theory of level spacing statistics 
for mixed systems back in 1984 0, there has been a lot of confusion in the literature 
advocating various phenomenological models due to incompatibility of experimental 
or numerical data with the Berry- Robnik (BR) statistics (see a recent comment 0). 
BR distribution is built on a simple and clean assumption of a statistically inde- 
pendent superposition of partial subspectra consisting of regular or chaotic levels 
(following an old Percivals' idea || of classifying the quantum eigenstates of mixed 
systems as regular or chaotic) . The sequence of regular levels, associated to eigen- 
states whose phase space distribution functions (e.g. Wigner or Husimi) localize on 
regions of regular motion, is assumed to have Poissonian statistics, whereas the se- 
quences of chaotic levels, associated with eigenstates whose phase space distribution 
functions extend over chaotic components of classical phase space, are assumed to 
have GOE (or GUE if antiunitary symmetry is absent) statistics of ensembles of 
Gaussian random matrices. Further, it is crucial to note that the gap distribution 
E(S), the probability that unfolded energy interval of length S contains no levels, 
factorizes upon independent superposition of level sequences, so the 2-component 
BR distribution for a system with a single classically chaotic component of relative 
measure p2 and regular components of complementary measure p± = 1 — p2 reads 

Ebr(S) = E Poisson (p 1 S)E GOE (p 2 S). (1) 

Note that -Epoisson^) = exp(— S) while for Eqq E (S) no closed-form expression exists 
(for the exact infinitely-dimensional GOE), and we have to rely on various expan- 
sions (we recommend Pade approximation published in ||]). The more common 
nearest neighbour level spacing distribution P(S) is directly related to the gap dis- 
tribution, simply as Pbr(S) = d 2 Ebr(S) / 'dS 2 . 

However, for the validity of semiclassical BR formula, two conditions have to be 
satisfied, (i) The regular and irregular levels should not be correlated, i.e. the cor- 
responding (Wigner or Husimi) phase-space distributions should not overlap. This 
is true if the quantum resolution scale in phase space, h 1 ^ 2 (where % is the effective 
Planck's constant), is small enough to resolve the essential features of the structure 
of classical phase space: h 1 ^ 2 < (sizes of the main regular islands, widths of chaotic 
strips penetrating through regular islands, etc), (ii) The quantum relaxation time, 
i.e. the Heisenberg (break) time tbreak = 2nh/AE (where AE is the mean level 
spacing) should be larger than the classical ergodic time t erg on the chaotic com- 
ponent, tbreak > ^erg- When this is not true, one expects dynamical localization of 



2 



eigenstates inside the chaotic component Jy, [H], [LI], [12 . 

Note that the BR statistics are incompatible with level repulsion, namely -Prr(O) = 
1 — p\ ^ 0. If either (i) or (ii) is violated, one recovers level repulsion P(S —>■ 0) —>■ 0. 
Indeed, numerous numerical studies (pi [7], |T3[ and references therein) give phe- 
nomenological support to the fractional power-law level repulsion which is usually 
very well globally captured by the phenomenological Brody distribution |14| 



p B (S) = (/3+l)&^exp(-^ +1 ), b=[T(l + (P + l)- 1 )f +1 (2) 

In fact, even for a generic 2-dim toy system with a simple phase space structure 
(where (i) and (ii) have the largest chances to apply) having a small number of 
islands and well connected chaotic component, one may verify that (i) and (ii) are 
typically fulfilled only for sequential quantum numbers substantially larger than 
~ 10 6 - 10 7 |L0|. 



So it is not surprising that the 'ultimate semiclassical' BR statistics have so 
far been clearly demonstrated only in two toy systems: (1) in a rather abstract 
compactified standard map |TJ|, and (2) in a 2-dim semi-separable oscillator [JIB], |TU 



which is dynamically a generic system but geometrically somewhat special. Here 
we give the first clear numerical demonstration of BR statistics in a generic billiard 
system with a smooth boundary. We consider classical and quantum motion of a 
free particle moving inside a bounded planar region which has a shape of a smoothly 
deformed circle. Billiard domain is described by the following function r(0), giving 
the radial distance from the origin to the boundary as a function of the polar angle 

0, 

r(0) = 1 + acos(40). (3) 

For the purpose of this letter we choose the following value of deformation param- 
eter, a = 0.04, for which the classical phase space (plotted in a Poincare-Birkhoff 
coordinates on a boundary-section in figure 1) has regular regions with the total 
relative Liouville measure (not the area on SOS Jl6|) pf = 0.115 ±0.005. Note that 
numerical computation of measures of regular and chaotic components of phase 
space in mixed (KAM) systems converges very slowly with increasing discretization 
of the phase space [T7|], hence it is difficult to further reduce the error estimate 
5pf « 0.005. 

High-lying quantum eigenenergies, eigenvalues of the Schrodinger equation (V 2 + 
A; 2 )\l/jt(r*) = with Dirichlet b.c. on the boundary r = r(0), have been computed 
by means of extremely efficient scaling technique proposed by Vergini and Saraceno 
fL8f : Eigenstates are expanded in a basis of circular scaling functions (see also 
I2|, as opposed to plane waves used in original approach flIB|) 



M 

*k{r) = E a iMkr) sin(4/0). (4) 

Note that the billiard has been desymmetrized and here we consider only fully 
antisymmetric states with respect to 8-fold symmetry group of the billiard. The 
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coefficients ai are determined by minimizing a special positive quadratic form de- 
fined along the boundary of the billiard |TJ| . The dimension of the problem M = 



[(1 + a)k/A] + M evan escent is nearly optimal where few ten, typically M evanescent ~ 40, 
evanescent modes have been added in order to ensure convergence and accuracy of 
the computed energy levels. We should note that the scaling method of quantiza- 
tion of billiards is by far superior to other relevant methods, e.g. boundary integral 



method |T9| or Heller's plane wave decomposition pD |, since it yields a constant 
fraction (5% — 10%) of M oc k of accurate levels, with no risk of missing any, by 
solving a single generalized eigenvalue problem of dimension M. 

In figure 2 we show cumulative nearest neighbour level spacing distribution 
W(S) = Jo dsP(s) = (d/dS)E(S)-(d/dS)E(0) for the unfolded spectral stretches 
{e n = k 2 /32 + (1/8 + l/ir)k n ; k min < k n < k max } (for small a) containing about 6000 
consecutive levels each. In fact, we have computed several spectral stretches, the first 
in the near semiclassical regime 399.7 < k < 600.1 (containing 6220 levels), and the 
last in the far semiclassical regime 15999.707 < k < 16004.865 (containing 5168 lev- 
els) where the sequential quantum number is N w k 2 /32 + (1/8 + l/ir)k « 0.8 • 10 7 . 
Only for the last spectral stretch in the far semiclassical regime (k « 16000) we 
found statistically significant agreement with BR distribution (figures 2,3) where the 
quantal (best-fitting) parameter p± agrees very well with its classical value, namely 
Pi = 0.119. However, for smaller sequential quantum numbers, when we approach 
the near-semiclassical regime, we find substantial deviation from BR statistics and 
recover fractional-power law level repulsion f2T| , namely for the lowest spectral 
stretch (figures 2,3) at k ~ 500 we find almost statistically significant agreement with 
Brody distribution @ with exponent (3 = 0.46. Of course, the fit to BR distribution 
in the near semiclassical regime k « 500 and the fit to Brody distribution in the far 
semiclassical regime k « 16000 turned out to be highly statistically non- significant. 

In figure 3 we show deviations of numerical spacing distributions from the semi- 
classical BR distribution (for parameter p\ = 0.119 ~ pf) in fine detail, using a 
smooth U-transformation pifl of the cumulative level spacing distribution U(W(S)) — 
U(W B r(S)), where U(W) = (2/tt) arccos^l - W, against W(S). This statistical 
representation has a uniform expected statistical error SU(W) = l/(7T V / AiV) (where 
AiV is the number of levels in a spectral stretch) and a constant density of numerical 
points along the abscissa. One can see very clearly that in both cases, far and near 
semiclassical, the numerical distributions are fluctuating around theoretical BR and 
Brody distributions, respectively, within expected statistical error. 

Finally we wish to characterize long-range spectral correlations as well, so we 
consider the number variance £ 2 (£) = (N 2 ) L — (N) 2 L , i.e. the variance of the number 
of unfolded levels e n in an interval of length L. Since this is a linear statistic it should 
be additive upon statistically independent superposition of spectral subsequences 



2~2~|. According to assumptions (i) and (ii) one immediately arrives to the ultimate 



semiclassical formula for the number variance |2! 

S 2 (L) = S2 oisson ( Pl L) + S^ OE (p 2 L) (5) 
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where £p oisson (L) = L is the number variance of Poissonian level sequence, and 
Sqoe(-^) ~ (2/tt 2 ) ln(27rL) is the number variance of the spectrum of infinitely di- 
mensional GOE random matrix which is supposed to model chaotic levels. In figure 
4 we show £ 2 (£) for four spectral stretches, namely for k « 500, k « 2000, k w 8000, 
and k « 16000, and only the last in the far semiclassical regime agrees well with the 
formula @ (for parameter p\ = pf = 0.115) up to L = L* « 50. 



In this letter we have clearly demonstrated the validity of BR level spacing distri- 
bution in a generic smooth plane billiard system with mixed classical phase space, 
namely the quartic billiard. However, for insufficiently small semi-classical parame- 
ter h ~ N~ l l 2 , we demonstrated the existence of fractional-power law level repulsion 
which is (for sufficiently small energy ranges) globally very well captured by the phe- 
nomenological Brody distribution. Unfortunately, this is the regime which can only 
be observed in most experimental situations due to extremely high energy region of 
crossover to BR statistics. 

We should note that this particular KAM billiard system ((|3|) for a = 0.04) has 
quite simple phase space structure which is reflected in relatively low transition point 
(N ~ 10 ) to the ultimate semiclassical BR statistics. For example, in a well known 
quadratic or Robnik billiard, the phase space is much more complicated ]23| (smaller 
regular islands, more partial phase space bariers, cantori), and as a consequence, the 



transition to BR regime is shifted to much higher energies [24 
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Figure 1: Classical phase space for a = 0.04 in Poincare-Birkhoff coordinates: 
arc-length s and tangential (normalized) velocity v. We show a chaotic orbit with 
2000000 collisions with the boundary. 
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Figure 2: Cumulative nearest level spacing distribution W(S) for a stretch of 
5168 consecutive levels in the far semiclassical regime (k m 16000) (thick curve) 
and a stretch of 6220 consecutive levels in the near semiclassical regime (k rs 500) 
(thin curve). The first numerical curve is almost overlapping with theoretical best 
fitting BR distribution for p± = 0.119 (dashed curve), while the second numeri- 
cal curve agrees very well with the best fitting Brody distribution with exponent 
(3 = 0.46 (dot-dashed curve). For comparison we give Poisson and GOE integrated 
level spacing distributions (dotted curves). In the inset we plot the same data in 
the T-function representation [^TJ, T(S) = ln(— ln(l — W(S))) against lnS 1 , which 
transforms the Brody distributions (an hence also Poissonian and Wigner) to straight 
lines, and enhances the region of small spacings. 
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Figure 3: Fine detail deviations from Berry- Robnik distribution (for pi = 0.119) in 
a uniform U-function transformation plf : we plot U(W(S)) — U(Wbr(S)) against 
W(S). In the far semiclassical regime k « 16000 (5168 consecutive levels), the 
difference of U-functions (thick curve) lies within a band of expected statistical error 
SU (dashed lines), while in the near semiclassical regime k ph 500 (6220 consecutive 
levels), the difference of U-functions (thin curve) agrees very well with the difference 
of U-functions for the best fitting Brody distribution with exponent (3 = 0.46 (dash- 
dotted curve). 
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Figure 4: Number variance S 2 (L) for the four spectral stretches: for k pa 16000 
(thickest curve, 5168 levels), k pa 8000 (next thickest curve, 17300 levels), k pa 2000 
(second thinnest curve, 5100 levels), and k pa 500 (thinnest curve, 6220 levels). 
Dashed curve is the semiclassical formula (H) which indeed reproduces the far semi- 
classical numerical data (thickest full curve) quite well, for L < L* pa 50. For 
comparison we give the Poissonian and GOE curves (dotted). 
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